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ABSTRACT 

We present a general procedure to solve numerically the general relativistic magnetohydrodynamics 
(GRMHD) equations within the framework of the 3 + 1 formali sm. The work repor ted here extends 
our previous investigation in general relativistic hydrodynamics (|Banvuls et alJl997|l where magnetic 
fields were not considered. The GRMHD equations are written in conservative form to exploit their 
hyperbolic character in the solution procedure. All theoretical ingredients necessary to build up high- 
resolution shock-capturing schemes based on the solution of local Riemann problems (i.e. Godunov- 
type schemes) are described. In particular, we use a renormalized set of regular eigenvectors of the 
flux Jacobians of the relativistic magnetohydrodynamics equations. In addition, the paper describes 
a procedure based on the equivalence principle of general relativity that allows the use of Riemann 
solvers designed for special relativistic magnetohydrodynamics in GRMHD. Our formulation and 
numerical methodology are assessed by performing various test simulations recently considered by 
different authors. These include magnetized shock tubes, spherical accretion onto a Schwarzschild 
black hole, equatorial accretion onto a Kerr black hole, and magnetized thick accretion disks around 
a black hole prone to the magnetorotational instability. 
Subject headings: MHD - relativity - methods: numerical 



1. INTRODUCTION 

In several astrophysical scenarios both magnetic and 
gravitational fields play an important role in determining 
the evolution of the matter. In these scenarios it is a com- 
mon fact the presence of compact objects such as neutron 
stars, most of which have intense magnetic fields of order 
10 12 - 10 13 G, or even larger at birth, ~ 10 14 - 10 15 G, 
as inferred from studies of anomalous X-ray pulsars and 
soft gamma-ray repeaters l|Kouveliotou et allll998|) . In 
some cases, i.e. in the so-called magnetars, the mag- 
netic fields can be so strong to affect t he internal struc- 
ture of the star ((Bocauet et alJll995(l . In a different 
context, the most promising mechanisms for producing 
relativistic jets like those observed in AGNs and micro- 
quasars, and the ones conjectured to explain gamma- 
ray bursts involve the hydromagnetic centrifugal accel- 
eration of material from an accretion disk, or the ex- 
traction of rot ational energy from the ergosphere of a 
Kerr black hole (IPenrosdl969tlBlandford fc ZnaiekH977t 
Bla ndford fc Pavnejll982TT In addition, the differential 
rotation of the magnetized plasma in the disk is respon- 
sible of the magnetorotational instability, which plays an 
impor tant role in transportin g angular momentum out- 
ward l)Balbus fc HawlevllT99"H) . 

If the gravitational field is strong enough, as in the 
vicinity of a compact object, the Newtonian description 
of gravity is only a rough approximation and general 
relativity becomes nec essary. In such a th eory, the so 
called 3+1 formalism ijArnowitt et al.lll962|) has proved 
particularly useful for numerical simulations involving 
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time-dependent computations of hydrodynamical flows 
in curved spacetimes, either static or dyn amic. The in- 
terested reader is addressed to lFontl <(2003D and references 
therein for an up-to-date overview of the different ap- 
proaches that have been introduced during the years for 
solving the general relativistic hydrodynamics equations. 

On the other hand, the inclusion of magnetic fields 
and the development of mathematical formulations of 
the magnetohydrodynamic (MHD) equations in a form 
suitable for efficient numerical implementations is still in 
an exploratory phase, although considerable progress has 
already been achieved in the last few years. 

Numerical studies in special relativistic magneto- 
hydrodynamics (SRMHD) have b een undertaken by 
a grow i ng n u mber of authors ij Komissarovl Il999t 
Balsa ral 120011: IKoldoba et all 12002 iDel Zanna et alJ 
l2Tcj3nLe ismann etHd"Il2"005j) . In particular. IKomissarmi 
1)1999(1 . iBalsaral |)20oI[k and IKoldoba et all (|2002f) 
developed independent upwind high-resolution shock- 
capturing (HRSC) schemes (also referred to as Godunov- 
type schemes), providing the characteristic information 
of the corresponding system of equations, which is the 
crucial b uilding block in su ch t ype of s c hemes . In 
addition, Ko missaro 3 1(1999]) and IBalsaral l(200l() pro- 
posed a comprehensive sample of tests to validate nu- 
merical MH D codes in spec ial relativity (SR). Re- 
cently, IDel Zanna et alJ 1(2003(1 have developed a third 
order shock-capturing central scheme for SRMHD which 
sidesteps the use of Riemann solvers in the solution pro- 
cedure (see, e.g. iTorol l(1997j) for general definitions on 
HRSC schemes). Simulations of the morphology and dy- 
namics of magnetized relativistic jets with Godunov-typ e 
schemes have been reported bv iLeismann et all {2005). 
In addition, the exact solution of the Riemann problem 
in SRMHD, for some particular orientation of the mag- 
netic field and the fluid velocity field, has been obtained 
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bv lRomero et al.l <|2005ft . 

Correspondingly, 3+1 repr esentations of GRMHD 
were first an al yzed b y iSloan fe Smarr _ 
lEvans fe Hawlevl (^1988^ . iZhand (119891). lYokosawal 

(1993 ), and, more re cently b v iKoide et all (11 998+ 
ipe Villiers fc H^Wl jxm&. IBanmga.rte fc Shaniro 

( 2003 ). IGammie et alJ 1 2003ft. iKomissarov |2005) , 
iDuez et alJ l|2005ft and iShibata fc Sekiguchil (12005ft . 



Most of the existing applications to date are i n the field 
of bla c k hole accretion and jet formation. In lYokosawal 
l)1993l 11995ft the transport of energy and angular 
momentum in magneto-hydrodynamical accretion onto 
a rotating black hole was studied adoptin g Wilson's for - 
mulation for the hydrodynamic equations (|Wilsonlll979ft . 
conveniently modified to account for the magnetic terms. 
The magnetic induction equation was solved using the 
constrained transport method of lEvans fc Hawlevl 
Later 



(1988). Later on, Koide and coworkers performed 
the first MHD simulat i ons o f jet f ormation in general 
relativity (jKoide et all 119981 12000ft in the context of 
the Blandford-Payne mechanism. These authors solved 
the MHD equations in the test-fluid approximation 
(in the background geometry of Schwarzschild/Kerr 
spacetimes) using a second-order finite difference cen- 
tral scheme with nonlinear dissipation. Employing 
the same nu merical approach IKoide et alJ l|2002ft and 
IKoide] l|2003ft studied the validity of the so-called MHD 
Penrose process to extract rotational energy from a 
Kerr black hole by simulating the evolution of a rarefied 
plasm a with a uniform magnetic field. iKomissarovl 
(2005) has also recently investigated this topic find- 
ing evidence in favour of the extraction of rotational 
energy of the black hole by the B landford-Znajek 
mechanism (Bland ford fc Znaiekl 119771) but against the 
development of strong relativ istic outflows or jet s. The 
long term solution found by IKomissarovl l)2005ft shows 
properties which are significantly different from those 
of the short initial (transient) phase studied by IKoide! 
(2003). An additional astrophysical application in the 
context of electromagnetic extraction of energy from 
a Kerr black hole is repre sented by the analysis of 
iMcKinnev fc Gammid ll2004ft. who have compare d the 
analytic prediction of felandfor d fc Znaiekl l|1977ft with 
time evolution calculations. Finall y, two different groups 
(|De Villiers fc Hawlevl l2003albt IGammie et aP \2W^i 
have started programs to investigate the time-varying 
behaviour of magnetized accretion flows onto Kerr 
black holes, with great emphasis on the issue of the 
development of the magnet orotational instability in 
thick accretion disks (see also l^ kosaw a fc Inuil (|2005)). 
While iDe Villiers fc Hawlevl (|20n8afLT ) adopt a noncon- 
ser vative ( ZEUS-li ke) sc heme, the approach followed 
by IGammie et al.l l|2003ft is based on a conservative 



HRS C scheme, na mely the so-called HLL scheme of 
Marten et aP l|1983ft . 

To the light of the existing literature on the subject it 
is clear that astrophysical applications of Godunov-type 
schemes in general rel ativistic MHD have o nly very re- 
cently been reported (IGammie et aP l2003t IKomissarovl 
I200a IDuez et all 12005ft . Our goal in this paper is to 
present the evolution equations for the magnetic field 
and for the fluid within the 3+1 formalism, formulated 
in a suitable way to apply Godunov-type schemes based 
on (approximate) Riemann solvers. Our numerical pro- 



cedure uses two original ingredients. On the one hand, 
the code incorporates a local coordinate transformation 
to Minkowskian coordinates, si milar to the one dev eloped 
for relativistic hydrodynamics in lPons eTaOCHI), prior 
to the computation of the numerical fluxes. In this way, 
Riemann solvers designed for SRMHD can be straight- 
forwardly used in GRMHD calculations. We note that 
IKomissarovl (12005ft applies the same approach, using a 
H RSC scheme based o n the SR Riemann solver described 
in IKomissarovl l|1999ft and adapted t o general re l ativit y 
following the procedure laid out in iPons et alJ (l998). 
We present here, however, a number of tests assessing the 
feasibility of the approach. As a second novel ingredient, 
we use a renormalized set of right and left eigenvectors of 
the flux vector Jacobians of the SRMHD equations which 
are regular and span a complete basis in any physical 
state, including degenerate states. 

The organization of the paper is as follows. We start by 
introducing the mathematical framework in §2, including 
the essentials of the 3+1 formalism, the description of 
the magnetic field, the induction equation and the con- 
servation equations of particle number, and stress-energy 
tensor in conservative form. A brief analysis of the hy- 
perbolic structure of the GRMHD system of equations 
is given in §3. The numerical procedure to solve the 
equations is described in §4. Finally, in §5 we present 
the results of some numerical tests and applications in 
order to assess our formulation and methodology. The 
summary of our work is given in §6. Throughout the 
paper Latin indices run from 1 to 3 and Greek indices 
from to 3. Four-vectors are indistinctly denoted using 
index notation or boldface letters, e.g. i// 1 , u. We adopt 
geometrized units by setting c = G = 1. 

2. MATHEMATICAL FRAMEWORK 

2.1. The Eulerian observer in the 3+1 formalism 

In the 3+1 formalism the line element of the spacetime 
can be written as 

ds 2 = -{a 2 - frP l )dt 2 + 2j3 i dx i dt + -/ ij dx i dx j , (1) 

where a (lapse function), /3 l (shift vector) and 7^ (spa- 
tial metric) are functions of the coordinates t, x % . A 
natural observer associated with the 3+1 splitting is the 
one with four velocity n perpendicular to the hypersur- 
faces of constant t at each event in the spacetime. This is 
the so-called Eulerian observer 3 . The contravariant and 
covariant components of n are given by 



»" = -(i, -n 

a 



(2) 



and 



»„ = (-a, 0,0,0), (3) 

respectively. In spacetimes containing matter an addi- 
tional natural observer is the one that follows the fluid 
during its motion, also called the comoving observer, with 
four-velocity u. With the standard definition, the three- 
velocity of the fluid as measured by the Eulerian observer 
can be expressed as 



hlu^ 



u • n 



(4) 



• j In the Kerr metric this Eulerian observer is indeed the ob- 
server with zero azimuthal angular momentum (ZAMO) as mea- 
sured from infinity. 
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where — u ■ n = W is the relative Lorentz factor between 
u and n, while h^ v — g^ + n^n^ is the the projector onto 
the hypersuface orthogonal to n, whose spatial terms are 
given by hij = 7y . From Eq. 0} it follows that 



(5) 



while 



Ui/W. Note that the Lorentz factor sat- 



isfies the relation W = 1/ 1/(1 — v 2 ) = aw* , where 
v 2 = j i jV l v : ' is the squared modulus of the three- velocity 
of the fluid with respect to the Eulerian observer. 

2.2. Magnetic field evolution 

A complete description of the electromagnetic field in 
general relativity is provided by the Faraday electromag- 
netic tensor field F^ v . This tensor is related to the 
electric and magnetic field, and B* 1 , measured by 
a generic observer with four- velocity £/ M , as follows, 



^imvXS being the volume element, 



n 



t 



:[/^A(5], 



(0 



(7) 



where g is the determinant of the 4- metric (g = det g^) 
and [/i^A(S] is the completely antisymmetric Levi-Civita 
symbol. Both, E and B are orthogonal to U, E ■ U = 
B • U = 0. The dual of the electromagnetic tensor *F^ V 
is defined as 

*F» V = lv^ XS Fx S , (8) 



and in terms of the electric and magnetic field measured 
by the observer U is given by 

^ vX5 U\E s . 



(9) 



From these equations, E and B can be expressed in terms 
of the electromagnetic tensor and the four-velocity U as 
follows 



E^ = F^ V U V , 
B» = *F» V U V . 



(10) 
(11) 



In terms of the electromagnetic tensor, Maxwell's equa- 
tions are written as follows, 



V U *F> ±V = Q, 



(12) 
(13) 



where Vy stands for the covariant derivative and J 11 is 
the electric four-current. According to Ohm's law, the 
latter can be in general expressed as 



J<* = Pq u» 



oF^i 



(14) 



where p q is the proper charge density measured by the 
comoving observer and a is the electric conductivity. 
Maxwell's equations can be further simplified if one as- 
sumes that the fluid is a perfect conductor. In this case 
the fluid has infinite conductivity and, in order to keep 
the current finite, the term proportional to the conduc- 
tion current, F^ v u v , must vanish, which means that the 
electric field measured by the comoving observer is zero. 
This case corresponds to the so-called ideal MHD con- 
dition. We can take advantage of this condition to ex- 
press the electric field measured by the observer U as a 



function of the magnetic field B measured by the same 
observer and of the four- velocities and u^. Straight- 
forward calculations give 



E» = ^ vXS u v U x B s . 



(15) 



If we choose U as the four-velocity of the Eulerian ob- 
server, U = n, Eq. I|15|l provides 



E a = 0, 
E l = -arj 



Qijk 



VjB k , 



(16) 
(17) 



or, in terms of three-vectors, E = —v x B, where the 
arrow means that the vector lies in the 'absolute space' 
and the cross product is defined using the induced volume 



element in the absolute space r\ 



arj 



Oijk 



Using the 



above relations, the dual of the electromagnetic field can 
be written in terms of the magnetic field only 



u»B v - u v B» 



II" 



(18) 



and Maxwell's equations V*.F Ml/ = reduce to the 
divergence-free condition plus the induction equation for 
the evolution of the magnetic field 



di^yB* 



dx l 



= 0, 



/T)B j 



x B 



(19) 

(20) 

(21) 
(22) 



-(av> -/3 J )B 1 }}, 
or, in terms of three- vectors, 
V ■ B = 

^K^)=Vx[(--^) 

2.3. Conservation Equations 

Once we have established the magnetic field evolution 
equation in the ideal MHD case, we need to obtain the 
evolution equations for the matter fields. These equa- 
tions can be expressed as the local conservation laws of 
baryon number and energy-momentum. For the baryon 
number we have 

V„ J u = 0, (23) 

J being the rest-mass current, J M = pu^, where p denotes 
the rest-mass density. The conservation of the energy- 
momentum is given by 



V V T^ = 0, 



(24) 



where T^ u is the energy-momentum tensor. For a fluid 
endowed with a magnetic field, this tensor is obtained by 
adding the energy- momentum tensor of the fluid to that 
of the electromagnetic field: 

(25) 



Fluid 



EM 



When the fluid is assumed to be perfect, TtY is given 

by 



1 Fluid 



phu»u v + pg^ , 



(26) 



where g^ v is the metric, p is the pressure, and h is the 
specific enthalpy, defined by h = 1 + e + p/p, e being the 
specific internal energy. The fluid is further assumed to 
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be in local thermodynamic equilibrium, and there exists 
an equation of state of the form p — p(p, e) which relates 
the pressure with p and e. On the other hand, the energy- 
momentum tensor of the electromagnetic field can 
be obtained from the electromagnetic tensor, F, as fol- 
lows 

T£u = ^ (f>*F" x - lg^F xs F xs ^j • (27) 

Furthermore, from Eq. JBJ and exploiting the ideal MHD 
condition, the electromagnetic tensor can be expressed in 
terms of the magnetic field 6 M measured by the comoving 
observer as 

pnu = -rj^s Uxb ^ (28) 



and Eq. I|27[l can be rewritten as 

1 

2 



(29) 



where b 2 = b v b v and where the magnetic field four vector 
has been redefined by dividing it by the factor y/An. As 
a result, the total energy- momentum tensor, fluid plus 
electromagnetic field, is given by 



= ph*u»u u + P *g" u - bn v 



(30) 



where we have introduced the definitions p* = p + b 2 /2 
and h* = h + b 2 /p. Note that if we consistently define 
e* = e+b 2 /(2p), the following relation, h* = 1+e* +p* / p, 
is fulfilled. 

In order to write the evolution equations (|23p. (|24[) in 

a conservation form suitable for numerical applications, 
let us define a basis adapted to the Eulerian observer, 

e (A)={n,^}, (31) 

where di are the coordinate vectors that are tangent to 
the hypersurface i=const, and, therefore, n-<9,; = 0. This 
allows us to define the following five 4- vectors T^(a) '■ 



V {A) ={T(e (A) ,0,J}, 



A = Q,.. 



(32) 



J (A) = t M e (A)>-; 

Hence the above system of equations (|23H . (|24|l can be 
written as 

V„X>fo = (33) 

where the five quantities S(m 011 the right-hand side -the 
sources-, are 



s (A) ={T^V M e (A)l/ ,0} 



(34) 



The covariant derivatives of the basis vectors, V /J e(A) ! y, 
are obtained in the usual manner as 



_ de Wv 



(35) 



where T vtl are the Christoffel symbols, and 

e(0)» = -a5 0v , e {k)v = g kv = (0k,lkj)- (36) 

In a similar way t o the pure hydrodynamics case 
(|Banvuls et all 11997(1 . if we now define the following 
quantities measured by an Eulerian observer, 

D = -J v n" = pW (37) 

Sj = -T(n, e(i) ) = ph*W 2 Vj - ab% (38) 

T = T(n,n) = ph*W 2 -p* -a 2 (b ) 2 - D (39) 



i.e. the rest-mass density, the momentum density of the 
magnetized fluid in the j-direction, and its total energy 
density (subtracting the rest-mass density in order to 
consistently recover the Newtonian limit), respectively, 
the system of GRMHD equations can be written explic- 
itly in conservative form. Together with the equation for 
the evolution of the magnetic field as measured by the 
Eulerian observer, Eq. the fundamental GRMHD 

system of equations can be written in the following gen- 
eral form 

1 /<VtF° dy/^gF^ 



-fJ 



dx° 



da 



(40) 



where the quantities F p (F° being the state vector and 
F 1 being the fluxes) are 



D 

Sj 

T 



Dv l 



S 3 v l +p*5) 

TV 1 + P *V l - 

v i B k - 



- bjB'/W 
v k B i 



(41) 



(42) 



with v l = v % — — . The corresponding sources S are given 

by 





Qk 



(43) 



where fe = (0,0, 0) T . Note that the following funda- 
mental relations hold between the four components of 
the magnetic field in the comoving frame, 6 M , and the 
three vector components B 1 measured by the Eulerian 
observer: 



WB l v, 
a 

B l + ab°u z 
W 



(44) 
(45) 



Finally, the modulus of the magnetic field can be written 

as 

32 



B 2 + a 2 ( b ) 2 
W 2 



b 2 = 



(46) 



where B 2 = B i B l . 



3. HYPERBOLIC STRUCTURE 



In Section |2~3*1 we have written the GRMHD equations 
in conservative form anticipating the use of numerical 
methods specifically designed to solve conservation equa- 
tions, as will be explained in the next Section. These 
methods strongly rely on the hyperbolic character of the 
eq uations and o n the associated wave structure. Follow- 
ing lAnilel |l989), in order to analyze the hyperbolicity of 
the equations it is convenient to write them in a more 
suitable form. If we take the following set of variables, 
V = (m m , b^jp, s), where s is the specific entropy, the sys- 
tem of equations can be written as a quasi-linear system 
of the form 

A<" B A V^V B = 0, (47) 
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where, A and B run from to 9, as the number of vari- 
ables, and the 10 x 10 matrices «4 M are given by 

-b»5% + P^bp l a » a ^\ 

0^ 



A" 



phS* 



\ 



U /3 







M 
0'' 







(48) 



where c s stands for the speed of sound 

. 

e being the mass-energy density of the fluid e = p(l ■ 
In Eq. (|48|l the following definitions are introduced: 

C = ph + b 2 , 

r"= 9 a " + 2!iY, 

1"° = {phg" a + (ph - b 2 /c 2 s )u>*u a )/ph 
f»°< = (u a b»/c 2 s -u IJ -b a )/ph, 
as well as the notation 

0^ = 0, 0" M ee (0,0,0, 0) T 



0^ ee (0,0,0,0). 



(49) 



(50) 
(51) 
(52) 
(53) 

(54) 



If (j){x^) = defines a characteristic hypersurface of the 
above system l|47[l. the characteristic matrix, given by 
A^de can be written as 



A e <b c = 



ph4> v 
V 0„ 



0„ a/ C 2 
Q„ a 



(55) 



where M = V p 0, a = u^(j) 



,„ B = W^, I" = = 
b 2 /(? s )au^/ph + Bb^/ph, /" = = 
V+2aui*)b u -B5Z. The 
determinant of the matrix l|55|l must vanish, i.e. 



^ + (ph 

(ab^/ci-Bu^/ph, and m£ = 



det(^ M ) = C a 2 A 2 JV 4 = 



(56) 



where 

A = Ca 2 

A/4 = ph 

and G = tfi^tfifj,. If we now consider a wave propagating 
in an arbitrary direction x with a speed A, the normal to 
the characteristic hypersurface is given by the four- vector 



-B, 
1 



a 2 G 



(57) 
B 2 G ,(58) 



0„ - (-A, 1,0,0), 



(59) 



and by substituting Eq. (|59")) in Eq. ffity we obtain the 
so called characteristic polynomial, whose zeroes give the 
characteristic speed of the waves propagating in the in- 
direction. Three different kinds of waves can be obtained 
according to which factor in equation l|56[) becomes zero. 
For entropic waves a = 0, for Alfven waves „4 = 0, and 
for magnetosonic waves A/4 = 0. 

Let us next analyze in more detail the characteristic 
equation. First of all, since the four-vector <^ must be 
spacel ike (this is a property of the RMHD system of equa- 
tions (|Anilei ll989)). it follows that ^(j)^ > 0. In terms 
of the wave speed A we obtain 



The characteristic speed A of the entropic waves prop- 
agating in the x-direction, given by the solution of the 
equation a = 0, is the following 



A = av x — (3 X 



(61) 



For Alfven waves, given by A — 0, there are two solu- 
tions corresponding, in general, to different speeds of the 
waves, 

b x ± VCu x , s 

A= 62 

b° ± VCu* K ' 

In the case of magnetosonic waves it is however not pos- 
sible, in general, to obtain explicit expressions for their 
speeds since they are given by the solutions of the quar- 
tic equation A/4 = with a, B and G explicitly written 
in terms of A as 



W 



a=—(-X + av x - (3 X ), 
a 

B = b x -b°X, 
1 



G = -(-(A + /3 1 ) 2 + tt 2 f 1 )- 



(63) 
(64) 
(65) 



(3 X < A < qa 



(60) 



Let us note that in the previous discussion about the 
roots of the characteristic polynomial we have omitted 
the fact that the entropy waves as well as the Alfven 
waves appear as double roots. These superfluous eigen- 
values appear associated with unphysical waves and are 
the result of working with the unco nstrained, 1 x 10 
system of equations. We note that Ivan Puttenl l)1991f) 
derived a different augmented system of RMHD equa- 
tions in constrained-free form with different nonphysi- 
cal waves. Any attempt to develop a numerical proce- 
dure based on the wave structure of the RMHD equa- 
tions must remove these nonphysical waves (and the 
corresponding eigenvectors) fro m the wave decom posi- 
tion. In the case of SRMHD Komissarov ( 1999J) and 
IKoldoba et alJ l|2002|) eliminate the nonphysical eigen- 
vectors by demanding the waves to preserve the values 
of the invarian ts u M u M = — 1 and u^b,,, = as su ggested 
bv lAnild <fl989L Correspondingly. iBalsaral 1)200 IS) selects 
the physical eigenvectors by comparing with the equiva- 
lent expressions in the nonrelativistic limit. 

It is worth noticing that just as in the classical case, 
the relativistic MHD equations have degenerate states in 
which two or more wavespeeds coinci de, which bre a ks the 
strict hyperbolicity of the system. iKomissarovl l)1999f) 
has reviewed the properties of these degeneracies. In 
the fluid rest frame, the degeneracies in both classical 
and relativistic MHD are the same: either the slow and 
Alfven waves have the same speed as the entropy wave 
when propagating perpendicularly to the magnetic field 
(Degeneracy I), or the slow or the fast wave (or both) 
have the same speed as the Alfven wave when prop- 
agating in a dir ection aligne d with the magnetic field 
(Degeneracy II) . I Anton et alJ l)2005t) have characterized 
these degeneracies in terms of the components of the 
magnetic field four-vector normal and tangential to the 
Alfven wavefront, b„, b t . When b„ = 0, the system 
falls within Degeneracy I, while Degeneracy II is reached 
when b t = 0. Let us note that the previous characteri- 
zation is covariant (i.e. defined in terms of four-vectors) 
and he nce can be checked in any reference frame. In ad- 
dition, lAnton et al.l l)2005f) have also worked out a single 
set of right and left eigenvectors which are regular and 
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span a complete basis in any physical state, including 
degenerate states. The renormalization procedure can 
be understoo d as a relativ i stic g eneralization of the work 
performed bv lBrio fc Wul |l988) in classical MHD. This 
procedure avoids the ambiguity inherent to a change of 
basis when appro aching; a degeneracy, as done e.g. by 
iKomissarovl (^999). The renormalized eigenvectors have 
been used in all the tests reported in the present paper 
using the full-wave decomposition Riemann solver. 

4. NUMERICAL APPROACH 

Writing the GRMHD equations as a first-order, flux- 
conservative, hyperbolic system allows us to use numer- 
ical methods specifically designed to solve such kind of 
equations. Among these methods, high-resolution shock- 
capturing (HRSC) schemes are recognized as the most 
efficient schemes to evolve complex flows accurately, cap- 
turing the discontinuities which appear when dealing 
with nonlinear hyperbolic equations. 
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4.1. Integral form of the GRMHD equations 

To apply HRSC techniques to the present GRMHD 
system we use Eq. (|40|) in integral form. Let f2 be a 
simply connected region of the four-dimensional mani- 
fold bounded by a closed three-dimensional surface dQ,. 
We take 30 as the standard-oriented hyperparallelepiped 
made up of the two spacelike surfaces S t , Sj+a* plus 
timelike surfaces £ x i , S^^ax* j that connect the two tem- 
poral slices. Then, the integral form of Eq. l|4U|l is 



1 5V7F° 



1 3^gF l 



iqV-9 9x° " ' Jn^f^g 3x l 
which can be written, for numerical purposes, as follows 



dVL= I SdO, 

n 

(66) 
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~^F 1 dx°dx 2 dx 3 
~gF 2 dx°dx 1 dx 3 
gF 3 dx°dx 1 dx 2 
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s<m, 



(67) 



where 



AV Li 



i'+Ai 1 r x 2 +Ax 2 r x 3 +Ax 3 



and 



y /jF°dx 1 dx 2 dx 3 
(68) 



AV = 



-Ax z 



y /jdx 1 dx 2 dx 3 



(69) 

The carets appearing on the fluxes denote that these 
fluxes, which are calculated at cell interfaces where the 
flow conditions can be discontinuous, are obtained by 
solving Riemann problems between the corresponding 
numerical cells. These numerical fluxes are further dis- 
cussed in Section 1431 

We note that in order to increase the spatial accuracy 
of the numerical solution, the primitive variables (see 
Sect. I4.5fl are reconstructed at the cell interfaces before 
the actual computation of the numerical fluxes. We use a 
standard second order minmod reconstruction procedure 
to compute the values of p, p, Vi and B % (i = 1,2,3) at 
both sides of each numerical interface. However, when 
computing the numerical fluxes along a certain direction, 
we do not allow for discontinuities in the magnetic field 
component along that direction. Furthermore, the equa- 
tions in integral form are advanced in time using the 
method of lines in conjunction with a second order, con- 
servative Runge-Kutta method l)Shu fc Osherlll988|) . 



4.2. Induction equation 

The main advantage of the above numerical procedure, 
Eq. (|67|l . to advance in time the system of equations, is 
that those variables which obey a conservation law are, 
by construction, conserved during the evolution as long 
as the balance between the fluxes at the boundaries of 
the computational domain and the source terms are zero. 
This is an important property that any hydrodynamics 
code should fulfill. 

However, as far as the magnetic field components are 
concerned, the system of equations (|4U|) only includes 
the induction equation Eq. i|22|). expressed by lj4*U|) 
in conservation form, while the divergence-free condi- 
tion, Eq. (|19|) . remains as an additional constraint to 
be imposed. Therefore, the numerical advantage of us- 
ing Eq. (|67|l for the conserved variables does not apply 
straightforwardly for the magnetic field components. In- 
deed, there is no guarantee that the divergence is con- 
served numerically when updating the magnetic field if 
we were to use the same numerical procedure we employ 
for the rest of components of the state vector. 

Among the methods designed to preserve the diver- 
gence of the magnetic field we use the constr a ined t rans- 
port method designed by IE vans fc Hawlevl 119881) and 
first ext ended to HRSC metho d s bv IRvu et alJ (flQQSfl 
(see also lLondrillo fc del Zannal l)2004|) for a recent dis- 
cussion) . This scheme is based on the use of Stokes the- 
orem after the integration of the induction equation on 
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surfaces of constant t and x l , St iX i. Let us write Eq. j^H 
as 



1 dB - ^ 

V7 9t 



(70) 



where we have defined the density vector B = \PyB and 

Q = (av-f3)x B. 

To obtain a discretized version of Eq. I|7()|l. we pro- 
ceed as follows. At a given time, each numerical cell 
is bounded by 6 two-surfaces. Consider, for concrete- 
ness, the two-surface T, t x 3, defined by t = const, and 
ir 3 = const., and the remaining two coordinates span- 
ning the intervals from x 1 to x 1 + Ax 1 , and from x 2 to 
x 2 + Ax 2 . The magnetic flux through this two-surface is 
given by 

(71) 



B ■ dY,. 



Furthermore, the electromotive force £ around the con- 
tour <9(£ t 3.3) is defined as 



S(t) = f 
Jd 



Qidx\ 



(72) 



Integrating Eq. (|70|l on the two-surface S t ^.3, and apply- 
ing Stokes theorem to the right hand side we obtain the 
equation 



dt 



-£ 



Vlidx 1 , 



which can be integrated to give 

r t+At 



(73) 



4> 



t+At 



fifda;* dt, (74) 



where the caret denotes again that quantities Cli are cal- 
culated at the edges of the numerical cells, where they 
can be discontinuous. At each edge, as we will describe 
below, these quantities are calculated using the solution 
of four Riemann problems between the corresponding 
faces whose intersection defines the edge. However, ir- 
respective of the expression we use for calculating flj, 
the method to advance the magnetic fluxes at the faces 
of the numerical cells satisfies, by construction, the di- 
vergence constraint. To see this we can integrate over a 
computational cell the divergence of the magnetic field at 
a given time. After applying Gauss theorem, we obtain 



V • BdV 



AV 



B-dY, 



6 

E 

faces, i— 



(75) 



In the previous expression, AV stands for the volume 
of a computational cell, whereas X denotes the closed 
surface bounding that cell. The summation is extended 
to the six faces (coordinate surfaces) shaping E. Now, 
taking the time derivative of Eq. I|75|l yields to 



— / \7-BdV = 
dt J AV 



E > 



dt 

faccs,i— 1 
6 4 

faces, i—1 edges, j — 1 



(76) 



where £ij is the contribution from edge j to the total 
electromotive force around the contour defined by the 
boundary of face i. It turns out that the above summa- 
tion cancels exactly since the value of £ for the common 
edge of two adjacent faces has a different sign for each 
face. Therefore, if the initial fluxes through each face of 
a numerical cell verify £?_ 



faces.' 



L $i = 0, this condition 
will be fulfilled during the evolution. 

4.3. Numerical fluxes and divergence-free condition 

The numerical integration of the GRMHD equations, 
Eqs. lj44lj) or JSJJ, is done using a HRSC scheme. Such 
schemes are specifically designed to solv e nonlinear hy- 
perbolic sy stems of conservation laws (|LeVecniel Il998t 
lTo7^IT997Ti . They are written in conservation form and 
use approximate or exact Riemann solvers to compute 
the numerical fluxes between neighbour grid zones. This 
fact guarantees the proper capturing of all discontinu- 
ities which may arise naturally in the solution space of 
a nonlinear hyperbolic system. Applications of HRSC 
schemes in rela ti vistic hydrodynami cs can be found in 
iMarti fe MulTerl <|2003|) : iFontJ {2003). Incidentally, we 
note that a detailed description of linearized Riemann 
S olvers based on the spectral decomposition can be found 
in lFont et alJ <I1994|) for special r elativistic hydrodynam- 
ics, and in Banvuls et alJ (119971) (diag onal metrics) and 
iFont et alJ l|2000D . Ilbanez et all l|2001li (general metrics) 
for general relativistic hydrodynamics. For HRSC meth- 
ods in classical MHD, on the other hand, we address to 



IRvu et alJl(T995llT99r^ . 

As discussed in Section |3 the existence of degenera- 
cies in the eigenvectors of the RMHD system of equa- 
tions makes it hazardous to implement linearized Rie- 
mann solvers based on the full spectral decomposition 
of the flux vector Jacobians. Nevertheless, we have suc- 
ceeded in developing and implementing in the code a full- 
wave decomposition (Roe-type) Riemann solver based on 
a single, renormalized set of right and left eigenvectors, 
as discussed in detail in lAnton et alJ (|2005T> . which is reg- 
ular for any physical state, including degeneracies. This 
Riemann solver is invoked in the code after a (local) lin- 
ear coordinat e transformation based on the procedure 
developed bv iPons et al.1 l)1998f) that allows to use spe- 
cial relativistic Riemann solvers in general relativity, and 
which has been properly extended to include magnetic 

fields (see Sect. 

In addition to the Roe-type Riemann solver we also 
use two simpler alternative approaches to compute the 
numerical flux es, namely the HL L single-state Rie- 
mann solver of lHarten et alJ <l!983j) and the second or- 
der ce ntral (symmetric) scheme of iKnrganov fc Tadmorl 
(2000) (KT hereafter). The KT scheme has proved 
recently to yield results with an accuracy compara- 
ble to those provided by full-wave decomposition Rie- 
mann solvers in simulations involvi ng purely hydrody- 
nami cal special relativistic flows l)Lucas-Serrano et alJ 
|2004|) and general re lativistic flows in dyn amical neu- 
tron star spacetimes llShibata fc Fontll 2005D. The inter - 
ested reader is addressed to .Kurganov fc TadmoH pOOOl: 
iLucas-Serrano et al.l l)2004|) for specific details on the KT 
central scheme. 

Correspondingly, the HLL Riemann solver is based on 
the calculation of the maximum and the minimum left 
and right propagating wave speeds emanating at the in- 
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terface between the two initial states, and the resulting 
flux is given by 

F(U L ,U fi ) = 

A + F(U L ) - A_F(U fl ) + A + A_(U^ - U £ ) 
A+ + A_ ' 

where A± = \±/a. Quantities F stand for the numeri- 
cal fluxes along each of the three spatial coordinate di- 
rections, namely F l (i — 1,2,3) in Eq. (|4(J|I . whereas 
U = F° denotes the state vector. Subindices L and R 
stand for the left and right states defining the Riemann 
problems at each numerical interface. Moreover A_ and 
A + are upper bounds of the speeds of the left- and right- 
propagating waves emanating from the cell interface, 

A + =max(0,A+ mSjL ,A+ mM? ), (78) 

A_=mm(0,A- msL ,A M ), (79) 

where Af ms I stands for the wavespeed of the fast mag- 
netosonic wave propagating to the left (s = —) or to the 
right (s = +) computed at state / (= L,R). These 
speeds are obtained by looking for the smallest and 
largest solution of the quartic equation A/4 = and can 
be effectively computed with a Newton-Raphson itera- 
tion scheme starting from A = ±0^/7" — P l (i = 1, 2, 3). 

Any of the flux formulae we have discussed can be 
used to advance the hydrodynamic variables according to 
Eq. (|67fl and also to calculate the quantities fti needed to 
advance in time the magnetic fluxes following Eq. (|74() . 
At each edge of the numerical cell, fii is written as an av- 
erage of the numerical fluxes calculated at the interfaces 
between the faces whose intersection define the edge. Let 
us consider, for illustrative purposes, Cl x . If the indices 
(J, k, I) denote the center of a numerical cell, an x— edge 
is defined by the indices (j, k+1/2, 1+1/2). By definition, 
Q x = a(v y B z - v z B y ). Since 

F y (B z ) — v y B z — v z B y (80) 

and 

F z (B y ) — v z B y — v y B z , (81) 
we can express £l x in terms of the fluxes as follows 

&xj,k+l/2,l+l/2 = 4 [-^fc+l/2,1 + Fj,k+l/2,l+l 

~ F j,k,l+l/2 ~ Fj,k+\,l+l/2h ( 82 ) 

where F y (F z ) refers to the numerical flux in the y (z) 
direction corresponding to the equation for B z (B y ) and 
multiplied by a to account for the correct definition of 
fl. Also note that in the numerical implementation of the 
constraint transport method, a sligh tly different proce- 
dure can be followed ijRvu et al.ll998tK According to this 
procedure, in the computation of the numerical fluxes 
<|80l) and (|8I|I . only the terms advecting the magnetic 
field are considered (i.e. the first term on the rhs of (|80[1 - 
(jHU), while the average in Eq. (|82|l is obtained dividing 
by a factor 2 instead of 4. Both of these procedures, the 
one described th rough Eqs. 18UI) - (|82|I and its modifica- 
tion provided bv IRvu et alJ ljl998j) allow us to advance 
the magnetic flux at the faces of the numerical cells in 
the correct way. However, we have also noted that for 



2D numerical tests our implementation of this modified 
scheme is gene rally more ro bust. We address the inter- 
est ed reader to iTotbl (|^000) for additional properties of 
the IRvu et aTl l)1998f) scheme. 

However, we need also to know the value of the mag- 
netic field at the center of the cells in order to obtain 
the primitive variables after each time step (cf. Sect. l4.5|l 
and to compute again the numerical fluxes of the other 
conserved variables for the next time step. If BJ ±1 , 2 k ; 
is the x-component of the magnetic field at the interface 
(j± 1/2, k, I), then the ^-component of the magnetic field 
at the center of the (j, k, I) cell, B^ k z , is obtained by tak- 
ing the arithmetic average of the corresponding fluxes, 
i.e. 

B j,k,l = jjC^ 3 'j-i/^,k,l^j-i/2,k,i + (83) 

where ASj ±1 ^ 2 k l is the area of the interface surface be- 
tween two adjacent cells, located at Xj±i/2 and bounded 
between [yk-i/2,yk+i/2] and [^-1/2,^+1/2]- Analogous 
expressions for ^ i+1 / 2 ,fc,m/2 and j+i/2,k+i/2,h and 
B y k l and B z k l can be easily derived. 

4.4. Special relativistic Riemann solvers in GRMHD 

In lPons et al.l l|1998j) we presented a general procedure 
to use any Riemann solver designed for the special rela- 
tivistic hydrodynamics equations in a general relativistic 
framework. In this section we describe a generalization 
of this approach to account for the magnetic field. It will 
be used to compute the numerical fluxes from the spe- 
cial relativistic full-wave decomposition Riemann solver 
discussed above. The procedure is based on performing 
linear transformations to locally flat (or geodesic) sys- 
tems of coordinates at each numerical cell interface, from 
which the metric becomes locally Minkowskian (plus sec- 
ond order terms). Notice that this approach is equiva- 
lent to the usual approach in classical fluid dynamics 
where one uses the solution of Riemann problems in slab 
symmetry for problems in cylindrical or spherical coor- 
dinates. 

In order to generalize this procedure to the GRMHD 
case one must start remembering that in the pure hy- 
drodynamical case, the components of the shift vector 
transversal to the cell interface play the role of a grid 
velocity, i.e., as if we have a moving interface. As dis- 
cussed in detail in IPons et al.l l|1998|h this can be eas- 
ily understood by noticing that, the fluxes through the 
moving interface for the local observer can be written as 

F % — h_ F°, where F % are the fluxes when 3 l = and F° 
the corresponding state vector. In terms of D, Sj, r and 
p* , the structure of the first five flux components (|42|l in 
the magnetic case follow the previous discussion with the 
conserved quantities advected with v l (that includes the 
correction term for the moving grid) and extra terms in 
the fluxes of momentum and energy (which do not de- 
pend explicitly on the shift vector ). This allows one to 
proceed along the same steps as in IPons et alJ l|1998() : i) 
Introduce the locally Minkowskian coordinate system at 
each interface; ii) solve the Riemann problem to obtain 
the numerical fluxes through the moving grid as seen by 
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the locally Minkowskian observer; iii) invert the trans- 
formation to obtain the numerical fluxes in the original 
coordinates. 

Let us now concentrate in the last three components of 
the fluxes Ij42(l . namely v l B k — v k B l , corresponding to the 
evolution of the magnetic field. The terms v % B k — v k B l 
also follow the discussion for the non-magnetic case and 
the same numerical procedure can be then applied. How- 
ever, the term f3 k B l /a couples the components of the 
shift vector parallel to the cell interface to the perpen- 
dicular magnetic field. This term has to be interpreted 
as a correction to the total electromotive force caused by 
the movement of the surface with respect to the Eulerian 
observer and has to be added to the final expression for 
the flux. 

In Section the validity of this approach with a full- 
wave decomposition Roe- type Riemann solver is assessed 
in a series of tests including discontinuous initial value 
problems, steady flows, and dynamical accretion disks. 
As a result of this assessment we conclude that the gen- 
eralized procedure to use SR Riemann Solvers in multidi- 
mensional GRMHD is an efficient and robust alternative 
to develop specific solvers that need of the knowledge of 
the whole spectral decomposition (eigenvalues and eigen- 
vectors) in general relativity. Since each local change of 
coordinates is linear and it only involves a few arithmeti- 
cal operations, the additional computational cost of the 
approach is negligible. 

4.5. Primitive variable recovering 

The numerical procedure used to solve the GRMHD 
equations allows us to obtain the values of the conserved 
variables F° at time t + At from their values at time t. 
However, the values of the physical variables (i.e. p, e, 
etc) are also needed at each time step in order to com- 
pute the fluxes. It is therefore necessary to solve the 
algebraic equations relating the conserved and the phys- 
ical variables. For the classical MHD equations and an 
ideal gas equation of state the physical variables can be 
expressed as explicit functions of the conserved ones. Un- 
fortunately, this cannot be done in GRMHD, a feature 
shared by the special and general relativistic versions of 
the purely hydrodynamics equations wit hin the 3+1 ap- 
proach (see iPapadopoulos fc Fontl (^99) for an alterna- 
tive formulation without this shortcoming). Therefore, 
the resulting nonlinear algebraic system of equations has 
to be solved numerically The procedure we describe be- 
low is an e xtension to full gen eral relativity of that de- 
veloped by iKomissarovl l)1999l) in the special relativistic 
case. 

The basic idea of this procedure relies on the fact that 
it is not necessary to solve the system H37|) - (Kffll) f° r the 
three components of the momentum, but instead for its 
modulus S 2 = S l Si. The next step is to eliminate the 
components of b a through Eqs. (|44Jl -H45 [l . After some 
algebra it is possible to write S 2 as 



(Z + B 2 ) 2 



W 2 



1 



(2Z 



(84) 



W 2 v ' Z 2 

where Z = phW 2 . 

The equation for the total energy can be worked out 
in a similar way 

B 2 (B l S^ 2 



Equations ({g4"|l and (JHSJ, together with the def- 

inition of Z, form a system for the unknowns p, p and 
W, assuming the function h — h(p,p) is provided. In 
our calculations we restrict ourselves to both, an ideal 
gas equation of state (EOS), p = pe(j — 1), for which 
h = 1 + 7p/p(7 — 1), where 7 is the adiabatic index, and 
a polytropic EOS (valid to describe isoentropic flows), 
p = Kp 7 , where K is the polytropic constant. In this 
last case the integration of the total energy equation can 
be avoided and the equation for the specific enthalpy is 
given by 



h = 1 



7' 



(86) 



Z + B^ 



p 



2W 2 



2Z 2 



D. 



(85) 



Then Eqs. P7)l and (|S4*Jl are solved to obtain p and W. 
5. RESULTS 

We turn now to assess the formulation of the GRMHD 
equations we have presented as well as the numerical 
techniques we employ to solve them. The simulations 
reported in this section are introduced in a way which 
gradually increases the level of complexity of the flow to 
solve, starting first with shock tube tests in both purely 
Minkowski spacetime and flat spacetimes suitably mod- 
ified by the presence of artificial gauge terms. Next we 
turn to one-dimensional tests of accreting magnetized 
flows onto Schwarzschild and Kerr black holes, to finally 
discuss two-dimensional simulations of thick accretion 
disks orbiting around black holes. This collection of tests 
allows us to validate our approach by comparing the nu- 
merical simulations with analytic solutions (in the cases 
where such a comparison is possible) , by investigating the 
ability of the code to preserve stationary solutions in the 
strong gravitational field regime, and by comparing with 
available numerical results reported in the literature. 

For those tests which involve (background) black hole 
spacetimes we adopt Boyer-Lindquist coordinates and we 
fix the unit of length to r g = M, M being the mass of 
the black hole. 

5.1. Relativistic Brio-Wu shock tube test 

The first test is the relativistic a nalog of the clas- 
sical Brio-Wu shock tube problem l)Brio & Wulli"983 
lBalsarall2001lh as ad apted to the relativistic MHD case 
bv lvan Putter] l|1993f l. The computational setup consists 
of two constant states which are initially at rest and sepa- 
rated through a discontinuity placed at the middle point 
of a unit length domain. The two states are characterized 
by the following initial conditions: Left state: p = 1.0, 
v x = 0.0, yy = 0.0, p = 1.0, and B« = 1.0. Right state: 
p = 0.125, v x = 0.0, yy = 0.0, p = 0.10, 5 y = -1.0. The 
adiabatic index of the ideal gas EOS is 7 = 2, and the 
x component of the magnetic field is equal for both left 
and right states, B x = 0.5. The test is performed using 
a Cartesian grid with 1600 cells. Results are reported for 
the HLL Riemann solver (as the other two schemes yield 
similar results) and for a CFL parameter equal to 0.5. 

The results of the simulation are shown in Fig.^ which 
displays the wave structure for various quantities after 
the removal of the membrane. This wave structure com- 
prises a fast rarefaction wave, a slow compound wave 
(both moving to the left), a contact discontinuity, and, 
moving to the right, a slow shock wave and a fast rar- 
efaction wave. The short dashed line in the six panels of 
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Fig. 1. — Wave pattern of the relativistic version of the Brio-Wu shocktube test. The left panel reports the rest-mass density p, the 
specific internal energy e and the y component of the magnetic field B v , while the right panel reports the x and y components of the 
velocity v x and v v , and the Lorentz factor W. The short dashed line shows the solution at time t = 0.4 in Minkowski spacctime. The 
open circles represent the solution at time t = 0.2 in Minkowski spacetime with gauge effects mimicked by a lapse function a = 2.0. The 
open squares represent the solution at time f = 0.4 in Minkowski spacetime with a shift vector /3 X = 0.4, while the long dashed line is the 
translation of the short dashed one by the amount (3 x t = 0.16. (Only 160 of the 1600 data points used in the simulation are drawn). 



Fig. ^ shows the wave pattern produced in the purely 
Minkowski spacetime at time t = 0.4. It is in good 
overal l agreement with the results obtained by iBalsaral 
(2001), in particular regarding the location of the differ- 
ent waves, the maximum value achieved by the Lorentz 
factor (W = 1.457), and the smearing of the numerical 
solution. In addition to this solution we use open cir- 
cles to denote the results of this test in flat spacetime 
but incorporating gauge effects by selecting a value of 
the lapse function different from unity, namely a = 2. 
The solution, which is shown at t = 0.2, matches as ex- 
pected with that represented by the short dashed line, 
obtained in flat spacetime at time t = 0.4. Finally, the 
open squares refer to a third version of this test carried 
out in a flat spacetime with a nonvanishing shift vector, 
namely (3 X = 0.4. The numerical displacement that is 
thus produced is in perfect agreement with the expected 
one. This is emphasized in the figure by translating the 
short dashed line into the long dashed one by the pre- 
dicted amount, (3 x t = 0.16. 

5.2. Magnetized spherical accretion 

In the second test we check the ability of the code to 
numerically maintain with a time-dependent system of 
equations the stationarity of the spherically symmetric 
accretion solution of a perfect fluid onto a Schwarzschild 
black hole in the presence of a radial magnetic field. It is 
worth emphasizing that a consistent solution for magne- 
tized spherical accretion with a force-free magnetic field 
satisfying the whole set of Maxwell equations does not 
exist (see Appendix lAlfor a proof). However, it is easy to 
show that any magnetic field of the type b a = (&', b r , 0, 0) 



does not affect the spherically symmetric hydrodynam- 
ical solution. Therefore, although the resulting config- 
uration is nonphysical, it provides a useful numerical 
test and has been used in the literature for this pur- 
pose Kl*mm\e et al.ll200li IT)e Villiers fc Hawlevll2003al 
IDuez et a,lJl2nn3T^ 

The initial setup consists of a perfect isoentropic fluid 
obeying a polytropic EOS with 7 = 4/3. The critical 
radius of the solution is located at r c — 8.0 and the rest 
mass density at the critical radius is p c — 6.25 x 10~ 2 . 
These parameters suffice to provide the full description 
of the spherical accretio n onto a Schw arzschild black hole 
as described in detail bv lMichell l|1972fl . The radial mag- 
netic field component, which can in principle follow any 
radial dependence, is chosen to satisfy the divergence- 
free condition. Moreover, its strength is characterized 
by the ratio (3 = b 2 /2p between the magnetic pressure 
and the gas pressure, computed at the critical radius of 
the flow. These initial conditions are evolved in time us- 
ing the Roe-type Riemann solver described in Sec. 14.41 
on a uniform radial grid covering the region between 
fmin = ^horizon + <5 and r max = 10.0, where 8 varies from 
0.1 to 0.3. 

Figure [21 shows the comparison between the analytic 
solution (solid lines) and the numerical solution (circles) 
for one representative case with pressure ratio = 1.0 
and 6 = 0.3. These results are obtained with a nu- 
merical grid of N — 100 radial zones, for which con- 
vergence is reached at time t = 250M. The order of 
accuracy of the code is computed by monitoring the er- 
ror L ee J2i=i \Qi - Qa,i\l YliLi Qa,i for quantity Q = p 
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Fig. 2. — Magnetized spherical accretion. Comparison between the analytic solution (solid line) and the numerical solution obtained 
with the Roe- type Riemann Solver (circles) for a model with j3 = 1.0 and N = 100, once convergence is reached. The left panels display 
the radial profiles of the rest mass density p (top) and the specific internal energy e, while the right panels show the corresponding profiles 
for the coordinate velocity v r (top) and the radial component of the magnetic field B r , all of them in geometrized units. 



as the number of grid points N is increased, where Q a 
represents the analytic solution. This procedure is re- 
peated for different values of the ratio /?, namely for 
(3 = 0, 1, 10, 100, and 1000 and the results, which are 
reported in Fig. [3] show that the global order of con- 
vergence of the code is 2, irrespective of the parameter 

A comparison of the accuracy of the three methods we 
use to compute the numerical fluxes is reported in Ta- 
ble^ f° r P = 10-0 and N = 70 radial zones. The results 
for the magnetized spherical accretion test appear in the 
upper half of the table. This table reports the global 
error of some representative quantities when numerical 
convergence is reached. For the particular test discussed 
in this section we find that there is not a single method 
providing the smallest error in all of the quantities, and 
the Roe-type Solver, which is the most accurate in the 
computation of the hydrodynamic variables, is the least 
accurate in the computation of the magnetic field. 

5.3. Equatorial Kerr accretion 

A further one-dimensional test of the code is provided 
by the stationary magnetized inflow solut ion in the Kerr 
metric derived bvlTakahas hi et ahl (^990). This solution 
was subsequently adapted to the case of equatorial in- 
flow in the region between the black ho l e horiz on and 
the marginally stabl e orbit bv |Gammid (1199 9) . This 
test has been u sed bylDe Villiers fc Hawlevl (2003a) and 
IGammie et"afl (|2003^ in the validation of their GRMHD 
codes. It represents a step forward in the level of com- 
plexity of the equations to solve with respect to those 
used in the previous two sections, since the test in- 
volves the Kerr metric, albeit specialized to the equa- 
torial plane. As a result, additional terms due to the in- 
creased number of nonvanishing Christoffel symbols ap- 




0.0001 



100 



1000 



N 



Fig. 3. — Error L of the rest mass density (see text for definition) 
for the magnetized spherical accretion test when the grid resolution 
is increased. The short-dashed and long-dashed lines indicate first 
and second order of global convergence, respectively. The sym- 
bols denote different values of the magnetization parameter at the 
critical point. 



pear in the equati ons. 

As described bv IGammid l|1999|) and adopting his no- 
tation, the inflow solution is determined once four con- 
served quantities are specified, namely the accretion rate 
Fm , the angular momentum flux Fl , the energy flux Fe 



12 



Q. 0.5 




0.3 



0.2 




3 

r/ r„ 



-0.006 



0.008 - 



-0.01 



0.012 - 



0.01 - 



0.008 




Fig. 4. — Comparison between the analytic solution (solid line) and the numerical converged solution obtained with the Roe- type 
Riemann Solver (circles) for the magnetized accretion solution onto a Kerr black hole with spin parameter a = 0.5. The left panel reports 
the rest mass density p and the azimuthal velocity tr , while the right panel reports the radial and the azimuthal components of the 
magnetic field, B r and B'r, all of them in geometrized units. 
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Fig. 5. — Error L of the rest mass density, the radial velocity and 
the toroidal magnetic field for the magnetized equatorial accretion 
in the Kerr metric. The dashed and the long-dashed line indicate 
first and second order of global convergence, respectively. 



and the component Fg$ of the electromagnetic tensor, 
which is related to the magnetic flux through the inner 
edge of the disk. For the sake of comparison we con- 
sid er an initial setup wit h the same numerical values used 
bv lGammie et all (|2003T> . namely a Kerr black hole with 
spin parameter a = 0.5, F M = -1.0, F L = -2.815344, 
F E = -0.908382, F e4> = 0.5. 



The numerical grid consists of N r x Ng gridpoints in 
the radial and angular directions, respectively. The ra- 
dial grid covers the region between r m - m = ^horizon + 0.2 
and r max = 4.0, while the angular grid consists of Ng = 3 
gridpoints subtending a small angle of 10 _5 7r accross the 
equatorial plane. The radial profiles of some significant 
variables, obtained with the Roe-type Riemann solver, 
are reported in Fig. 0] for a radial grid of N r = 100 
zones. The open circles indicate the numerical results 
while the underlying solid lines correspond to the ana- 
lytic solution. It is found that the stationarity of the 
solution is preserved to high accuracy by the numerical 
code. For the long-term evolutions considered there are 
no significant deviations from the analytic profiles. 

As we did for the magnetized spherical accretion test 
we use the current test to compute again the order of 
convergence of the code as the grid is refined. The global 
order of convergence for some representative quantities 
is reported in Fig.[5j which shows that the code is second 
order accurate. As already commented bv lGammie et all 
( 2003), the worsening of the order of convergence for 
at high grid resolution is due to the fact that the initial 
condition is "semi-analytic" , requiring the solution of an 
algebraic equation. Thus, the inaccuracies produced at 
time t = become more pronunced for large numbers of 
radial zones N r . 

The performance of the code using the HLL and KT 
solvers has also been checked with this test. While the 
order of convergence is preserved irrespective of the nu- 
merical shemes used to compute the fluxes, the actual 
accuracy can vary significantly. The results of this com- 
parison for the equatorial Kerr accretion solution are 
summarized in the lower half of Table ^ which reports 
the global error of representative quantities, when con- 
vergence is reached, on a numerical grid with N r = 60 
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TABLE 1 

Comparison among different schemes 



bp Sv r 8B r SB^ 

Michel test 

HLL 3.76 x 10~ 3 3.92 x 10~ 3 7.64 x lCT 1 ' = 

Roe-type 2.97 x 10" 3 3.45 x 10~ 3 1.09 x 10~ 12 

KT 3.36 x 1CT 3 3.54 x 10~ 3 1.94 x IP" 18 - 

Gammie test 

HLL 1.92 x 10~ 2 2.54 x 10" 3 2.28 x 10" y 1.48 x 1CT 3 

Roe-type 6.90 x 10" 3 3.01 x 10" 3 3.96 x 10" 3 2.14 x 10" 3 

KT 1.63 x 10" 2 9.72 x 10" 4 2.30 x 10" 9 9.89 x 10~ 3 



Note. — Accretion tests: Comparison of the accuracy of some repre- 
sentative quantities for the HLL, Roe and KT solvers. The columns report 
the global errors when convergence is reached. 



radial points. It is worth stressing that the HLL scheme, 
at least in our implementation, turns out to be the most 

I 



5.4. Thick accretion disks around black holes 

An intrinsic two-dimensional test for the code is 
provided by the stationary solution of a thick disk 
(or torus) orbiting a r ound a black hole, descr i bed b y 
iFishbone fc MoncrieJ (I1976T). IKozlowski et all l)1978jl . 
and more recently bv lFont &: Daiend l)2002j) . The result- 
ing configuration consists of a perfect barotropic fluid in 
circular non-Keplerian motion around a Schwarzschild or 
Kerr black hole, with pressure gradients in the vertical 
direction accounting for the disk thickness. These thick 
disks may posses a cusp on the equatorial plane through 
which matter can accrete onto the black hole. 

In the following two subsections we describe our nu- 
merical tests for unmagnetized and magnetized thick 
disks, respectively. In both cases the effective potential 
at the inner edge of the disk is smaller than that at the 
cusp, thus providing initial conditions which are strictly 
stationary. For simplicity we limit our simulations to 
models with constant distribution of specific angular mo- 
mentum t = —Uip/ut, although the same qualitative re- 
sults have been obtained with more general rotation laws. 

5.4.1. Unmagnetized disk 

In testing the evolution of a purely hydrodynamical 
torus we consider a model s imilar to the one used by 
iDe Villiers &: Hawlevl l)2003a(l for the Schwarzschild met- 
ric, namely a torus with specific angular momentum £ = 
4.5, position of the maximum density at r contor = 15.3, 
and an effective potential at the inner edge such that 
the inner and outer radii on the equatorial plane are 
r in = 9.34 and r out = 39.52, respectively. We choose a 
polytropic EOS with 7 = 4/3 and a polytropic constant 
K such that the torus-to-hole mass ratio is M t /M ~ 0.07. 

We have checked that the code can keep the stationar- 
ity of the initial equilibrium torus when evolved in time. 
Figure [S] shows the global order of convergence as com- 
puted from the rest mass density p. The correspond- 



accurate in the computation of the magnetic field. 



r 




0.001 L- 1 — 1 — '— 1 1 1 1 1 

100 

N 

Fig. 6. — Unmagnetized thick accretion disk. Error L of the 
rest mass density when resolution is increased. The short-dashed 
and the long-dashed lines indicate first and second order of global 
convergence, respectively. 



ing global error L reported in the figure, and defined as 

L = Z)fi=i \Pij ~ Pa,ij\/Y,i,j=i Pa,ij, is computed after 
10 orbital periods for each model, using a uniform numer- 
ical grid consisting of N x N gridpoints, whose specific 
values can be read off from the figure. As it is apparent 
from Fig.Elthe code reaches second order of convergence 
for reasonable high values of N (> 200). 

We note that in addition to the model just discussed 
we have also analyzed the performance of the code by 
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20 10 20 30 40 10 20 30 40 

Fig. 7. — Velocity field and logarithimic isocontours of the rest mass density. At four orbital periods an elongated high density structure 
is formed near the equatorial plane, signaling the MRI instability in 2D calculations. 



comparing the evolution of additional hydrodynam ical 
models which were studied bv lFont fc Daignd <|2002l) and 
iZanotti et al.l l)2003t) using independent codes based on 
HRSC schemes. In all the cases considered, correspond- 
ing to a number of different generalizations such as disks 
with power-law distributions of the specific angular mo- 
mentum, disks in Kerr spacetime, and disks subject to 
the so-called runaway instability, the GRMHD code re- 
produced the same quantitative results of the indepen- 
dent hydrodynamical codes with negligible differences. 

5.4.2. Magnetized disk 

As a final test we consider the evolution of a magne- 
tized torus around a Schwarzschild black hole. In this 
case, however, a stationary solution which might provide 
self-consistent initial data for such magnetized disks is 
not available. Indeed, it can be proved (see Appendix lAl 
for a proof) that the hydrodynamical isoentropic type 
of models that we have used in the previous section for 
unmagnetized disks cannot be "dressed" with a mag- 
netic field, to produce a force-free magnetized torus that 
satisfies the whole set of Maxwell's equations. There- 
fore, we follow the same prag mati c approach adop ted by 
iDe Villiers fe Hawlevl l|2003a|i and lGammie et alJ(|2003j) . 



and simply add an ad-hoc poloidal magnetic field to the 
hydrodynamical thick disk model. The magnetic field is 
generated by a vector potential oc max(p/p c — C, 0), 
where p c is the maximum rest mass density of the torus 
and C is a free parameter which determines the confine- 
ment of the field inside the torus. The hydrodynamical 
torus is the same as the one considered in Section l5.4.1l 
but endowed with a magnetic field characterized by a 
confinement parameter C — 0.5 and such that the aver- 
age ratio of magnetic-to-gas pressure inside the torus is 
=1.5xlCT 3 . 

The four panels of Fig.[3display isocontours of the rest 
mass density, logarithmically spaced, during the first few 
orbital periods of the evolution. These results correspond 
to a simulation employing the HLL solver with a compu- 
tational grid of 200 r adial zones and 10 angu lar zones. 
It was first shown by IBalbus fc Hawlevl <|1991^l that the 
dynamics of such magnetized thick disks is governed by 
the so-called magnetorotational instability (MRI) , which 
generates turbulence in the disk and helps explaining the 
transport of angular momentum outwards. In axisym- 
metry the development of the MRI is much less signif- 
icant than in full three dimensions and manifests itself 
through the appearence of the so-called "channel solu- 
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Fig. 8. — Top panel: the equatorial angular momentum ap- 
proaches the Keplerian profile as an effect of the magnetorotational 
instability. Bottom panel: time evolution of the total gas pressure 
(solid line) and of the magnetic pressure (dashed line). 



tion" <|De Villiers fc Hawlevll2003b|) . This feature of the 
solution becomes visible in our simulation after about 
three orbital periods, as shown in Fig. [7| in the form of 
a high density elongated structure near the equatorial 
plane. We report in Fig.[8]two additional distinctive fea- 
tures that can be unambiguosly attributed to the MRI. 
The first one, showed in the top panel, represents the 
transport of angular momentum (initially constant) out- 
ward, which acquires a Keplerian profile (indicated by a 
thick solid line) as the evolution proceeds. Correspond- 
ingly, the botton panel shows the rapid increase of the 
(mean) magnetic pressure (dashed line) with respect to 
the gas pressure (solid line) during the first two orbital 
periods and due to the MRI driven turbulence. 

We note, however, that the present status of the nu- 
merical code does not allow us to evolve efficiently ad- 
ditional simulations with higher resolutions and with in- 
creasingly larger values of the magnetization parameter. 
As a result, the typical distorsion of the isodensity con- 
tours produced by the MRI is not visible in Fig. [7| A 
parallel version of the code is currently under develop- 
ment. This will allow for higher resolution simulations 
of magnetized disks in astrophysical contexts. 

6. CONCLUSIONS 

In this paper we have presented a procedure to solve 
numerically the general relativistic magnetohydrody- 
namic equations within the framework of the 3 + 1 for- 
malism. The work reported here represents the extension 



of our previous investigation ijBanvuls et al.llT997j) where 
magnetic fields were not considered. The GRMHD equa- 
tions have been explicitely written in conservation form 
to exploit their hyperbolic character in the solution pro- 
cedure using Riemann solvers. Most of the theoretical 
ingredients which are necessary in order to build up high- 
resolution shock-capturing schemes based on the solution 
of local Riemann problems have been discussed. In par- 
ticular, we have described and implemented three alter- 
native HRSC schemes, either upwind as HLL and Roe, or 
symmetric as KT. Our implementation of the Roe-type 
Riemann solver has made use of the equivalence prin- 
ciple of general relativity which allows to use, locally, 
the characteristic information of the system of equations 
in the special relativistic limit, following a slight mod- 
ificatio n of the procedure first presented in iPons et alJ 
(1998). Further information regarding the renormaliza- 
tion of the eigenvectors of the GRMHD flux-vector Ja- 
cobians has been deferred to an accompanying paper 
l)Anton et alJ 12005). The work reported in this paper, 
hence, follows the recent stir of activity in the ongo- 
ing efforts of developing robust numerical codes for the 
GRMHD system of equations, as exemplified by the in- 
vestigations pre^ejrtedJii_tdT^jBst few years by a number 
of groups llDe Villiers fc Hawlevl l2003at iGammie et alJ 
200$ iDuez et al.ll2005t lKomissarovll2005D . 

Our formulation of the equations and numerical pro- 
cedure have been assessed by performing the various 
test simulations discussed in earlier works in the liter- 
ature, including magnetized shock tubes in flat space- 
times, spherical accretion onto a Schwarzshild black hole, 
equatorial magnetized accretion in the Kerr spacetime, as 
well as evolution of thick accretion disks subject to the 
development of the magnetorotational instability. The 
code has proved to be second order accurate and has 
successfully passed all considered tests. In the near fu- 
ture we plan to apply this code in a number of astro- 
physical scenarios involving compact objects where both 
strong gravitational fields and magnetic fields need be 
taken into account. 
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APPENDIX 

MAGNETIZED MICHEL ACCRETION 

In this Appendix we prove that there is not a consistent solution for a force-free magnetic field added to the 
spherically symmetric accretion of a perfect fluid onto a Schwarzschild black hole. In general, it is not at all obvious 
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that a hydrodynamical solution can be "dressed" with a force-free magnetic field. lOronl l)2002|) has shown that the 
form of the four-current compatible with a force-free magnetic field is given by 

J" = Pq U» + T)b» (Al) 

where p q is the proper charge density. Note that when r\ — 0, i.e. when the current is only due to the convective 
term, the assumption of force-free is automatically guaranteed by the ideal MHD condition. However, we will consider 
here the more general expression given by Eq. (|A1I) . If we write explicitely the four vanishing components of the 
electric field in the comoving frame of the accreting fluid, F^ v u v = 0, recalling that the velocity field is given by 
= (u° , u 1 , u 2 , u 3 ) = (u*,u r ,0, 0), we find 

F 01 =0, (A2) 
F 02 u° + F 12 u 1 = Q, (A3) 
F 3X =0, (A4) 

where we have also used the fact that F03 = <9o^3 — O^Aq = 0. Let us next consider the first couple of Maxwell 
equations 

F[ap,i\ = 0, (A5) 
where the comma denotes partial differentiation. After writing them explicitly for all possible combinations we obtain 

-Fbi,2 + -Fi2,o + i*20,i = 0, (A6) 
Fqi,3 + -Pi3,o + -^30,1 = 0, (A7) 

Fq2,3 + ^23,0 + -^30,2 = 0, (A8) 
^12,3 + ^23,1+^31,2=0. (A9) 

By the symmetries of the spacetime and by relations (|A2|) - (|A4|) this system reduces to 



^02,1=0, (A10) 
F 23 ,i=0. (All) 

Summarizing, among the 6 components of the antisymmetric electromagnetic tensor F^, 3 of them vanish, namely 
Fqi = Fq3 = F13 = 0. Among the rema ining 3, on ly two are independent, since the constraint l|A3|l has to be fulfilled. 
Furthermore, according to Eqs. l)A10jl and (|A11(I . F02 and F23 are functions of the angle only, F02 = Fq2(9) and 
F23 = ^23(^)1 and are therefore constants along fluid lines. Taking all this into account we can write the components 
of the magnetic field explicitly, using definition (|llfl in the main text 

6° = ^Lf 23 ui, (A12) 



-9 

^ = --L^o = ^ (A13) 

b 2 = Q, (A14) 

b 3 = —^=(Fq2Ui — F12M0) = 7=~T- ( Al 5) 

Note that Eq. (|A15|) can be alternatively computed from the condition b^u^ = 0. 

Up to this point we have shown that the magnetic field is completely determined by two constants, -F23 and -Fo2- 
We now consider the seco nd couple of Maxwell equations, namely V \ J F fJ - L/ — AnJ^. According to the assumption on 
the four-current, Eq. (|Alll . and on the four- velocity in the case of spherical accretion, these equations become 

d 2 (V^dF 02 )=47T^(p q u° + V b ), (A16) 

d 2 (^F 12 )=4n^d( Pq u 1 + V b 1 ), (A17) 

a 1 (7^F 21 ) = 0, (A18) 

d 2 (^/^gF 32 )^4T:^^r 1 b 3 , (A19) 

where F 02 = ^02/(700522 and F 12 — F\ 2 j ' g\\g 2 i- From l|A18|) it follows that the term Fy^lgw must be a function 
of the angular coordinate 9 only which, recalling (|A3|) and the fact that both u° and it 1 are functions of r, implies 
that F12 = Fq2 =0. As a result, the toroidal component of the magnetic field b 3 vanishes. Moreover, according to 
Eq. I|A19|) . the term i 7 2 3 /(^ 2 sin 9 ) mus t be a function of r only. Given that F23 = F2z{0), it must be F23 = Asin#, 
with A a constant. Finally, l|A16|l and i|A17|) are now reduced to the following homogeneous system in the unknowns 
p q and r\ 

u° Pq + b Q r] = 0, (A20) 
u x p q + fo 1 // = 0. (A21) 
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Imposing the vanishing of the determinant gives b° /b 1 — u° /u\ which cannot be satisfied since it violates the constraint 
coming from the combination of the orthogonality condition b^u^ — and the normalization condition u^u^ = — 1. 
This concludes the proof that it is not possible to add a force-free magnetic field to the hydrodynamic solution of 
spherical accretion in a Schwarzschild spacetime that satisfies the full set of Maxwelll equations. 



APPENDIX 

MAGNETIZED THICK ACCRETION DISK 

In this Appendix we show that it is not possible to build a consistent stationary and axisymmetric solution for a 
magnetized torus by simply adding a force-free magnetic field to t he hydrodynamic equilibrium model of an isoentropic 
thick accretion disk (Kozlowski et al.lll978i iFont fc Daignc 2003). The proof, that for simplicity we limit to the case 
of Schwarzschild spacetime but can be extended to a Kerr black hole as well, could follow the same reasoning of 
the previous Appendix. However, the demonstration is more direct if one exploits some topological properties of 
the expected solution. In fact, from Maxwell equations it is possible to show that the magnetic field of a perfectly 
conducting medium endowed with a purely toroidal motion has to be purely poloidal, i.e. b r ^ 0, b 6 ^ 0, while 
h 1 = b^ = 0. Under these conditions the magnetic field lines lie on the surfaces of constant magnetic potential 
(magnetic surfaces), which coincide with the surfaces of constant angular velocity O = u^/u*. This property prevents 
the generation of a toroidal component of the magnetic field, even in the presence of differential rotation (Ferraro's 
theorem), and allows to introdu ce a new coo rdinate system (x\,X2) such that x\ varies along the poloidal field lines 
and X2 is constant along them ijOronl 12002). In this new coordinate system the magnetic field will only have one 
non- vanishing component b 1 , while b 2 = . 

According to lBekenstein fc Oronl l)1979|) . for a force- free magnetic field in an isoentropic flow the quantity U = u l u^ 
is constant along the magnetic surfaces, and it can be used to define the new coordinate Xi ■ In the case of circular 
motion in Schwarzschild spacetime this quantity reads 



" = -SFTS) (A1 » 

where t — ~u^ju t is the specific angular momentum. According to von Zeipel's theorem ijvon Zeinell 9241) t is constant 
along surfaces of constant £1 for the class of barotropic hydrodynamic models that we are considering. Therefore, both 
and I are constant along magnetic surfaces and the new coordinate X2 can be defined as 

>f -(-*f-<i^. <-) 

which is the so-called von Zeipel parameter l)Chakrabartilll985j) . The other coordinate x\ can be chosen such that 
orthogonality between x\ and X2 is preserved, i.e. g\2 — 0. After some calculations involving straightforward metric 
coefficient transformations, this choice yields to 

Xy = (r-3M)cos0 . (A3) 

In computing Eq. (|A3|) we have made the reasonable ansatz that x% is factorized as X\ = p(r)q(8). lOronl ll200l has 
shown that, in order to satisfy the second couple of Maxwell's equations and the scalar equation V M (/i6 M ) = 0, which 
can be proved to hold for any isoentropic magnetized flow, the following factorization in terms of generic functions of 
X\ and X2 must exist 

siw =/(x,, ' !(I2) ' (A4) 

where A = —gtt9<p<t> m the Schwarzschild metric. From the normalization condition u^u^ = — 1 it follows that 
(u*) 2 = l/[(g tt (l - a^fl 2 )], and Eq. JUJl becomes 



2 

2 



(x 2 y(l - n\x 2 ) z )- z = f(x 1 )h(x 2 ). (A5) 



Since f2 = fl(x2), Eq. I|A5(1 requires that the term 1 — 2M/r is factorizable as f{x{)h{x2), which can be shown not to 
be possible. Hence, the constraint i|A4f) cannot be met, and a force-free magnetized torus built from the isoentropic 
hydrodynamic model of a thick accretion disk cannot be obtained. 
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